<!DOCTYPE html>
<html lang="en">
  <head>
    <meta charset="utf-8">
    <title>Cindy JS</title>
    <script type="text/javascript" src="../../build/js/Cindy.js"></script>
    <script type="text/javascript" src="../../build/js/CindyGL.js"></script>
    <link rel="stylesheet" href="../../css/cindy.css">
  </head>
    
  <body style="font-family:Arial;">
    
    <h1>CindyJS: Raycasting Surfaces</h1>
    
    <script id="csmousedown" type="text/x-cindyscript">
    sx = mouse().x;
    sy = mouse().y;
    dragging = sx < .5;
    </script>
    <script id="csmouseup" type="text/x-cindyscript">
    dragging = false;
    </script>
    <script id="csdraw" type="text/x-cindyscript">
    //the following is executed for every rendered frame
    if (dragging,
        dx = 3 * (sx - mouse().x); dy = 3 * (sy - mouse().y);,
        dx = .005 * cos(seconds() * .1); dy = .003 * sin(seconds() * .1);
    );
    
    sx = mouse().x;
    sy = mouse().y;
    
    //the rotation matrix: It is modified either if the user is dragging or time passes
    mat = mat * (
        (1, 0, 0),
        (0, cos(dy), -sin(dy)),
        (0, sin(dy), cos(dy))
    ) * (
        (cos(dx), 0, -sin(dx)),
        (0, 1, 0),
        (sin(dx), 0, cos(dx))
    );


    //the 3 sliders at the left.
    PA.x = 0.55;
    if (PA.y > .4, PA.y = .4);
    if (PA.y < -.4, PA.y = -.4);
    a = (.4 + PA.y) / .8 * 1.1;

    PB.x = 0.6;
    if (PB.y > .4, PB.y = .4);
    if (PB.y < -.4, PB.y = -.4);
    alpha = ((.4 + PB.y) / .8) * .7 + .3;

    PC.x = 0.65;
    if (PC.y > .4, PC.y = .4);
    if (PC.y < -.4, PC.y = -.4);
    zoom = exp(3 * PC.y - 1);
    
    //configuration for the lights in the scene. A light has a position, a gamma-parameter for its shininess and a color
    lightdirs = [
        ray((.0, .0), -100), //enlights parts of the surface which normal points away from the camera
        ray((.0, .0), -100),
        ray((.0, .0), 100), //Has an effect, if the normal of the surface points to the camera
        ray((.0, .0), 100),
        (-10, 10, -2.),
        (-10, 10, -2.),
        (10, -8, 3.),
        (10, -8, 3.)
    ];

    gamma = [2, 20, 2, 20, 1, 10, 1, 10];

    colors = [
        (.3, .5, 1.),
        (1, 2, 2) / 2,
        (1., 0.2, 0.1),
        (2, 2, 1) / 2,
        .4 * (.7, .8, .3),
        .9 * (.7, .8, .3),
        .4 * (.6, .1, .6),
        .9 * (.6, .1, .6)
    ];


    colorplot(computeColor(#)); //render the scene. # is the pixel coordinate
    
    
    drawtext((-.65, -.45), "degree: $" + if(newN<100,newN,"\infty") +"$");

    //lines for the sliders
    draw((.55, .4), (.55, -.4), color -> (0, 0, 0));
    draw((.6, .4), (.6, -.4), color -> (0, 0, 0));
    draw((.65, .4), (.65, -.4), color -> (0, 0, 0));    
    </script>
               
    <script id="csinit" type="text/x-cindyscript">
    //initialize some variables
    mat = [
        [0.3513, -0.4908, -0.7973],
        [-0.8171, -0.5765, -0.0051],
        [-0.4571, 0.6533, -0.6036]
    ];
    sx = mouse().x;
    sy = mouse().y;
    dragging = false;
    N = 5;
    zoom = 2.2;
    a = 0.3;
    alpha = .7;
    eps = .001;
    
    //we stand at position mat*(0, 0, -2.2/zoom) and watch to (0,0,0).
    //ray(pixel, t) is the point in R^3 that lies at position t the ray behind the pixel at location pixel(vec2)
    //t=0 is corresponds to points within the interesting area near (0,0,0)
    ray(pixel, t) := mat * ((t+2.2/zoom) * (pixel.x, pixel.y, 1) + (0, 0, -2.2/zoom));
    
    //sphere with radius 1/zoom for clipping
    S(r) := (r * r - 1/(zoom*zoom));
    
    //fun is the user defined trivariate polynomial
    fun(x, y, z) := (x ^ 2 + y ^ 2 + z ^ 2 - (0.5 + a) ^ 2) ^ 2 - (3.0 * ((0.5 + a) ^ 2) - 1.0) / (3.0 - ((0.5 + a) ^ 2)) * (1 - z - sqrt(2) * x) * (1 - z + sqrt(2) * x) * (1 + z + sqrt(2) * y) * (1 + z - sqrt(2) * y);
    
    //F takes vec3 instead of 3 variables
    F(p) := fun(p.x, p.y, p.z);
  

    //guess the degree of the trivariate polynomial F. This approximation is reliable up to degree ~20.
    helpfun(s, x) := log(|F(s*x)|)/log(s*|x|); //is approx. degree+log(leadingcoeff)/log(s*|x|) for large s
    guessdeg() := max(apply(1 .. 2, //take the best result of 2
      x = [random(), random(), random()];
      s = 1;
      l = 1;
      best = 1;
      it = 0;
      while(l<100 & s < 1e50 & it<100, //throw away Infinity
        best = round(l);
        it = it+1;
        s = 2*s;
        l = 2*helpfun(s*s, x)-helpfun(s,x); //remove error caused by log(leadingcoeff)
      );
      if(it==100, best = 1000000);
      best        
    ));
    
    
    //casteljau algorithm to evaluate and subdivide polynomials in Bernstein form.
    //poly is a vector containing the coefficients, i.e. p(x) = sum(0..N, i, poly_(i+1) * b_(i,N)(x)) where b_(i,N)(x) = choose(N, i)*x^i*(1-x)^(N-1)
    casteljau(poly, x) := (
      regional(alpha, beta);
      alpha = 1-x;
      beta = x;
      forall(0..N, k,
        repeat(N-k,
          poly_# = alpha*poly_# + beta*poly_(#+1);
        );
      );
      poly //the bernstein-coefficients of the polynomial in the interval [x,1]
    );

    //evaluates a polynomial, represented as vector of coefficients in bernstein-form
    eval(poly, x) := casteljau(poly, x)_1;
    
    //this function has to be called whenever N changes
    init() := (
        newN = guessdeg(); errc(newN);
        if(newN!=N,
          N = newN;
          //The following line is DIRTY, but it makes the application run smooth for high degrees. :-)
          //Nethertheless, it might cause render errors for high degree surfaces. In fact, only a subset of the surface is rendered.
          //Adapt limit according to hardware.
          //values of kind 4*n-1 are good values, as it means to use vectors of length 4*n.
          N = min(N,7); 
          
          //N+1 Chebyshev nodes for interval (0, 1)
          li = apply(1..(N+1), k, (cos((2 * k - 1) / (2 * (N+1)) * pi)+1)/2);
          
          //A is the matrix of the linear map that evaluates a polynomial in bernstein-form at the Chebyshev nodes
          A = apply(li, node,
            //the i-th column contains the values of the (i,N) bernstein polynomial evaluated at the Chebyshev nodes
            apply(0..N, i, eval(
              apply(0..N, if(#==i,1,0)), // e_i = [0,0,0,1,0,0]
              node //evaluate  b_(i,N)(node)
            )) 
          );
          
          B = (inverse(A)); //B interpolates polynomials (in Bernstein basis), given the values [p(li_1), p(li_2), ...]
        )
        
    );
    init();
    
    //B3 is a matrix that interpolates quadratic polynomials (in monomial basis), given the values [p(-2), p(0), p(2)]
    B3 = inverse(apply([-2, 0, 2], c, apply(0 .. 2, i, c ^ i))); 

    //use central difference to approximate dF
    dF(p) := (
        (F(p + [eps, 0, 0]) - F(p - [eps, 0, 0])),
        (F(p + [0, eps, 0]) - F(p - [0, eps, 0])),
        (F(p + [0, 0, eps]) - F(p - [0, 0, eps]))
    ) / (2 * eps);


    

    //update the color color for the pixel at position pixel assuming that the surface has been intersected at ray(pixel, dst)
    //because of the alpha-transparency updatecolor should be called for the intersections with large dst first
    updatecolor(pixel, dst, color) := (
      regional(x, normal);
      color = (1 - alpha) * color;
      x = ray(pixel, dst); //the intersection point in R^3
      normal = dF(x);
      normal = normal / abs(normal);

      forall(1..length(lightdirs),
        //illuminate if the normal and lightdir point in the same direction
        illumination = max(0, (lightdirs_# / abs(lightdirs_#)) * normal);
        color = color + alpha * (illumination ^ gamma_#) * colors_#;
      );
      color
    );
    

    
    nsign(pixel, a, b) := ( //Descartes rule of sign for the interval (a,b)
      //obtain the coefficients in bernstein basis of F along the ray in interval (a,b) by interpolation within this interval
      poly = B * apply(li,
        F(ray(pixel, a+#*(b-a))) //evaluate F(ray(pixel, ·)) along Chebyshev nodes for (a,b)
      );
      //count the number of sign changes
      ans = 0;
      //last = poly_1;
      forall(2..(N+1), k,
        //if(last == 0, last = poly_k;); this (almost) never happens
        if(min(poly_(k-1), poly_k) <= 0 & 0 <= max(poly_(k-1), poly_k), //sign switch; avoid products due numerics
          ans = ans + 1;
        );
      );
      ans //return value   
    );
    
    
    //bisect F(ray(pixel, ·)) in [x0, x1] assuming that F(ray(pixel, x0)) and F(ray(pixel, x1)) have opposite signs
    bisectf(pixel, x0, x1) := (
        regional(v0, v1, m, vm);
        v0 = F(ray(pixel, x0));
        v1 = F(ray(pixel, x1));
        repeat(11,
            m = (x0 + x1) / 2; vm = F(ray(pixel, m));
            if (min(v0,vm) <= 0 & 0 <= max(v0, vm), //sgn(v0)!=sgn(vm); avoid products due numerics
                (x1 = m; v1 = vm;),
                (x0 = m; v0 = vm;)
            );
        );
        m //return value   
    );
    
    
    //id encodes a node in a binary tree using heap-indices
    //1 is root node and node v has children 2*v and 2*v+1
    
    //computes s=2^depth of a node id: Compute floor(log_2(id));
    //purpose: id corresponds interval [id-s,id+1-s]/s
    gets(id) := (
      s = 1;
      repeat(10,
        if(2*s<=id,
          s = 2*s;
        )
      );
      s //return value
    );
    
    //determines the next node in the binary tree that would be visited by a regular in DFS
    //if the children of id are not supposed to be visited
    //In interval logic: finds the biggest unvisited interval directly right of the interval of id.
    next(id) := (
      id = id+1;
      //now: remove zeros from right (in binary representation) while(id&1) id=id>>1;
      repeat(10,
        if(mod(id,2)==0, 
          id = floor(id/2);
        )
      );
      if(id==1, 0, id) //return value - id 0 means we stop our DFS
    );
    
    //what color should be given to pixel with pixel-coordinate pixel (vec2)
    computeColor(pixel) := (
      regional(a,b,u,l);
        
        spolyvalues = apply([-2, 0, 2], S(ray(pixel, #))); //evaluate S along ray
        spoly = B3 * spolyvalues;                          //interpolate to monomial basis

        D = (spoly_2 * spoly_2) - 4. * spoly_3 * spoly_1; //discriminant of spoly
        
        color = gray(.7); //the color, which will be returned

        if (D >= 0, //ray intersects ball
            l = (-spoly_2 - re(sqrt(D))) / (2 * spoly_3); //intersection entering the ball
            u = (-spoly_2 + re(sqrt(D))) / (2 * spoly_3); //intersection leaving the ball

            a = l;
            b = u;
            
            //traverse binary tree (DFS) using heap-indices
            //1 is root node and node v has children 2*v and 2*v+1
            id = 1; 
            //maximum number of steps
            repeat(N*6,
              //id=0 means we are done; do only a DFS-step if we are not finished yet
              if(id>0,
                s = gets(id); //s = floor(log_2(id))
                
                //the intervals [a,b] are chossen such that (id in binary notation)
                //id = 1   => [a,b]=[l,u]
                //id = 10  => [a,b]=[l,(u+l)/2]
                //id = 101 => [a,b]=[l,(u+3*l)/4]
                //id = 11  => [a,b]=[(u+l)/2,u]
                //...
                a = u - (u-l)*((id+1)/s-1);
                b = u - (u-l)*((id+0)/s-1);
                
                //how many sign changes has F(ray(pixel, ·)) in (a,b)?
                cnt = nsign(pixel, a, b);
                if(cnt == 1 % (b-a)<.01*zoom, //in this case we found a root (or it is likely to have a multiple root)
                  //=>colorize and break DFS
                  color = updatecolor(pixel, bisectf(pixel, a, b), color);
                  id = next(id),
                if(cnt == 0, //there is no root
                  id = next(id), //break DFS
                  
                  //otherwise cnt>=2: there are cnt - 2*k roots.
                  id = 2*id;  //visit first child within DFS
                )
            );  
          ));
            
        );
        color //return value
    );
    </script>
    

    <div  id="CSCanvas" style="border:0px"></div>
    
    <script type="text/javascript">
        var cdy = CindyJS({canvasname:"CSCanvas",
                    scripts: "cs*",
                    animation: {autoplay: true},
                    use : ["CindyGL","katex"],
                    geometry: [ { name:"PA", kind:"P", type:"Free", pos: [.5,.3,1], narrow: true, color: [1,1,1], size:8 },
                                { name:"PB", kind:"P", type:"Free", pos: [.5,.2,1], narrow: true, color: [1,1,1], size:8 },
                                { name:"PC", kind:"P", type:"Free", pos: [.5,.1,1], narrow: true, color: [1,1,1], size:8 } ],
                    ports: [{
                      id: "CSCanvas",
                      width: 700,
                      height: 500,
                      transform: [ { visibleRect: [ -0.7, -0.5, 0.7, 0.5 ] } ]
                    }]
        });
    </script>
    
    <div>
<input type="text" id="inp" value="(x^2+y^2+z^2-(0.5+a)^2)^2-(3*((0.5+a)^2)-1)/(3-((0.5+a)^2))*(1-z-sqrt(2)*x)*(1-z+sqrt(2)*x)*(1+z+sqrt(2)*y)*(1+z-sqrt(2)*y)"  onkeypress="if((event.which ? event.which : event.keyCode)==13) { cdy.evokeCS('fun(x,y,z) := (' + this.value + '); init();'); }" size="80" style="font-size:18px">
<select id="sel" onchange="document.getElementById('inp').value = this.value; cdy.evokeCS('fun(x,y,z) := (' + this.value + '); init();');" style="width:20em;">  
  <option value="(x^2+y^2+z^2-(0.5+a)^2)^2-(3*((0.5+a)^2)-1)/(3-((0.5+a)^2))*(1-z-sqrt(2)*x)*(1-z+sqrt(2)*x)*(1+z+sqrt(2)*y)*(1+z-sqrt(2)*y)">Kummer Quartic</option>
  <option value="4*((a*(1+sqrt(5))/2)^2*x^2-y^2)*((a*(1+sqrt(5))/2)^2*y^2-z^2)*((a*(1+sqrt(5))/2)^2*z^2-x^2)-1*(1+2*(a*(1+sqrt(5))/2))*(x^2+y^2+z^2-1)^2">Barth Sextic</option>
  <option value="-2*a/125+x^8+y^8+z^8-2*x^6-2*y^6-2*z^6+1.25*x^4+1.25*y^4+1.25*z^4-0.25*x^2-0.25*y^2-0.25*z^2+0.03125">Chmutov Octic</option>
  <option value="a*(-1/4*(1-sqrt(2))*(x^2+y^2)^2+(x^2+y^2)*((1-1/sqrt(2))*z^2+1/8*(2-7*sqrt(2)))-z^4+(0.5+sqrt(2))*z^2-1/16*(1-12*sqrt(2)))^2-(cos(0*pi/4)*x+sin(0*pi/4)*y-1)*(cos(pi/4)*x+sin(pi/4)*y-1)*(cos(2*pi/4)*x+sin(2*pi/4)*y-1)*(cos(3*pi/4)*x+sin(3*pi/4)*y-1)*(cos(4*pi/4)*x+sin(4*pi/4)*y-1)*(cos(5*pi/4)*x+sin(5*pi/4)*y-1)*(cos(6*pi/4)*x+sin(6*pi/4)*y-1)*(cos(7*pi/4)*x+sin(7*pi/4)*y-1)">
    Endraß Octic
  </option>
  <option value="x^2+y^2+z^2-1">Ball</option>
  <option value="k = 6; x^k+y^k+z^k-1">Cube</option>
  <option value="x^2+z^2-1/3*(1+y)^3*(1-y)^3">Citric</option>
  <option value="x^2+y^2+z^3-z^2">Ding Dong</option>
  <option value="x^3+x^2*z^2-y^2">Hummingbird</option>
  <option value="x^2-x^3+y^2+y^4+z^3-z^4">Vis a Vis</option>
  <option value="(x^2+9/4*y^2+z^2-1)^3-x^2*z^3-9/80*y^2*z^3">Sweet</option>
  <option value="k=a*2;(x+(k/2-1))*(x^2+y^2+z^2-k^2/4)+z^2">Parabolic Ring Cyclide</option>
  <option value="cos(x)*sin(y) + cos(y)*sin(z) + cos(z)*sin(x)+a">Gyroid</option>
  <option value="cos(2*x)+cos(2*y)+cos(2*z)+a">Schwarz P</option>
</select>
</div>
    Non-algebraic functions are approximated by polynomials.
    Roots are isolated by Descartes Method in Bernstein basis.
  </body>
</html>
